The Astrophysical Journal 

Preprint typeset using BTpX style emulateapj v. 04/20/08 



LA-UR-09-06131 



o 

(N 
< 

o 

u 

43 
6 



(N 
> 
O 

Os 
O 



13 



THE COYOTE UNIVERSE III: SIMULATION SUITE AND PRECISION EMULATOR FOR THE NONLINEAR MATTER 

POWER SPECTRUM 

Earl Lawrence 1 , Katrin Heitmann 2 , Martin White 3 , David Higdon 1 , Christian Wagner 4 5 , Salman Habib 6 , and 

Brian Williams 1 

1 CCS-6, CCS Division, Los Alamos National Laboratory, Los Alamos, NM 87545 

2 ISR-1, ISR Division, Los Alamos National Laboratory, Los Alamos, NM 87545 

3 Departments of Physics and Astronomy, University of California, Berkeley, CA 94720 

4 Astrophysikalisches Institut Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam 

4 Institut de Ciencies del Cosmos (ICC), Universitat de Barcelona, Marti i Franques 1, E-08028 Barcelona and 
6 T-2, Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545 
Submitted to the Astrophysical Journal 

ABSTRACT 

Many of the most exciting questions in astrophysics and cosmology, including the majority of observational 
probes of dark energy, rely on an understanding of the nonlinear regime of structure formation. In order to fully 
exploit the information available from this regime and to extract cosmological constraints, accurate theoretical 
predictions are needed. Currently such predictions can only be obtained from costly, precision numerical 
simulations. This paper is the third in a series aimed at constructing an accurate calibration of the nonlinear 
mass power spectrum on Mpc scales for a wide range of currently viable cosmological models, including dark 
energy. The first two papers addressed the numerical challenges, and the scheme by which an interpolator was 
built from a carefully chosen set of cosmological models. In this paper we introduce the "Coyote Universe" 
simulation suite which comprises nearly 1,000 A^-body simulations at different force and mass resolutions, 
spanning 38 wCDM cosmologies. This large simulation suite enables us to construct a prediction scheme, 
or emulator, for the nonlinear matter power spectrum accurate at the percent level out to k ~ 1 /zMpc -1 . We 
describe the construction of the emulator, explain the tests performed to ensure its accuracy, and discuss how 
the central ideas may be extended to a wider range of cosmological models and applications. A power spectrum 
emulator code is released publicly as part of this paper. 

Subject headings: methods: A^-body simulations — cosmology: large-scale structure of universe 



1. introduction 

During the last three decades, cosmology has made tremen- 
dous progress: from order of magnitude estimates to measure- 
ments of key cosmological parameters approaching percent 
level accuracy. The standard model of cosmology, based on 
the growth of structure by gravitational instability, has been 
impressively validated on large scales where the perturbations 
to the Friedman model are small. While the model is success- 
ful at predicting or reproducing a wide array of observations, 
it contains several mysterious elements with perhaps the most 
mysterious being t he accelerated expan sion of the Universe 
dRiess et alj|l998t iPerlmutter et al.lll999l) . The accelerated 
expansion may be caused by a dark energy or hinting at a 
modification of general relativity on the largest scales. 

To date most of the best known cosmological param- 
eters have been constrained primarily from the study of 
anisotropies in the cosmic microwave background (CMB) ra- 
diation. However, large-scale structure probes play an im- 
portant role in breaking parameter degeneracies and in con- 
straining conditions in the late-time Universe. Such probes 
are becoming ever more precise, with future surveys aiming 
at measurements approaching percent level accuracy to bet- 
ter characterize the Universe in which we live. Techniques 
based on the observation and analysis of cosmic structure in- 
clude the use of baryon acoustic oscillations, redshift space 
distortions, weak lensing measurements and the abundance 
of clusters of galaxies; they stand to play a pivotal role in 
improving our understanding of the dynamics of the Uni- 
verse. (For a recent discussion regarding improvements on 



dar k energy constraints f r om co mbining different probes, see, 
e.g. lAlbrecht et alJ l2006. 2009). The large scale structure of 
the Universe contains information about both the geometry, as 
well as the dynamics of structure formation. In combination 
these two pieces of information can help distinguish between 
dark energy or a modification of general relativity as the prime 
cause of cosmic acceleration. 

On small scales, the cosmological interpretation of struc- 
ture formation probes is complicated due to the nonlinear 
physics involved. Co mmonly-used fitting funct i ons for e.g., 
the p ower spectrum dPeacock & DoddsTll996t jSmith et al. 
120031) have poorly characterized systematics and are no longer 
adequate for precision work. Absent a controlled theoret- 
ical framework, direct use of simulations (augmented with 
phenomenological parameters as appropriate) is essential if 
the physics is to be more correctly captured. The simula- 
tion codes need to be adequately tested to ensure they meet 
the new demands being placed upon them. The simulations 
which pass these goals are often very expensive and only a 
restricted number of runs can be performed. This in turn puts 
a premium on developing very efficient strategies to constrain 
parameters from limited observations and simulations. High- 
precision prediction schemes for different statistics are essen- 
tial to succeed in this task. 

This paper is the third in a series aimed at addressing this 
question in the context of the nonlinear matter power spec- 
trum on Mpc scales. Current observations in weak lensing 
are quickly becoming theory limited due to the lack of pre- 
cise theoretical estimates of this quantity for a wide range of 
cosmological models. This is a pressing problem. It is also 
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TABLE 1 

The parameters for the 37+1 models which define the sample space; iLi is measured in Mpc -1 . See text for further details. 
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M000 
M001 
M002 
M003 
M004 
MOOS 
M006 
M007 
M008 
M009 
MO 10 
M011 
M012 
M013 
MOM 
M015 
M016 
M017 
M018 



0.1296 
0.1539 
0.1460 
0.1324 
0.1381 
0.1358 
0.1516 
0.1268 
0.1448 
0.1392 
0.1403 
0.1437 
0.1223 
0.1482 
0.1471 
0.1415 
0.1245 
0.1426 
0.1313 



0.0224 
0.0231 
0.0227 
0.0235 
0.0227 
0.0216 
0.0229 
0.0223 
0.0223 
0.0234 
0.0218 
0.0234 
0.0225 
0.0221 
0.0233 
0.0230 
0.0218 
0.0215 
0.0216 



0.9700 
0.9468 
0.8952 
0.9984 
0.9339 
0.9726 
0.9145 
0.9210 
0.9855 
0.9790 
0.8565 
0.8823 
1.0048 
0.9597 
1.0306 
1.0177 
0.9403 
0.9274 
0.8887 



1.000 
0.816 
0.758 
0.874 
1.087 
1.242 
1.223 
0.700 
1.203 
0.739 
0.990 
1.126 
0.971 
0.855 
1.010 
1.281 
1.145 
0.893 
1.029 



0.8000 
0.8161 
0.8548 
0.8484 
0.7000 
0.8226 
0.6705 
0.7474 
0.8090 
0.6692 
0.7556 
0.7276 
0.6271 
0.6508 
0.7075 
0.7692 
0.7437 
0.6865 
0.6440 



0.7200 
0.5977 
0.5970 
0.6763 
0.7204 
0.7669 
0.7040 
0.6189 
0.7218 
0.6127 
0.6695 
0.7177 
0.7396 
0.6107 
0.6688 
0.7737 
0.7929 
0.6305 
0.7136 



0.12 
0.11 
0.10 
0.11 
0.14 
0.12 
0.14 
0.11 
0.12 
0.13 
0.12 
0.13 
0.15 
0.14 
0.14 
0.14 
0.13 
0.13 
0.14 



0.19 
0.18 
0.17 
0.17 
0.22 
0.20 
0.24 
0.18 
0.20 
0.21 
0.19 
0.21 
0.24 
0.23 
0.23 
0.22 
0.20 
0.21 
0.23 



MO 19 
M020 
M021 
M022 
M023 
M024 
M025 
M026 
M027 
M028 
M029 
M030 
M031 
M032 
M033 
M034 
M035 
M036 
M037 



0.1279 
0.1290 
0.1335 
0.1505 
0.1211 
0.1302 
0.1494 
0.1347 
0.1369 
0.1527 
0.1256 
0.1234 
0.1550 
0.1200 
0.1399 
0.1497 
0.1485 
0.1216 
0.1495 



0.0232 
0.0220 
0.0221 
0.0225 
0.0220 
0.0226 
0.0217 
0.0232 
0.0224 
0.0222 
0.0228 
0.0230 
0.0219 
0.0229 
0.0225 
0.0227 
0.0221 
0.0233 
0.0228 



0.8629 
1.0242 
1.0371 
1.0500 
0.9016 
0.9532 
1.0113 
0.9081 
0.8500 
0.8694 
1.0435 
0.8758 
0.9919 
0.9661 
1.0407 
0.9239 
0.9604 
0.9387 
1.0233 



1.184 
0.797 
1.165 
1.107 
1.261 
1.300 
0.719 
0.952 
0.836 
0.932 
0.913 
0.777 
1.068 
1.048 
1.147 
1.000 
0.853 
0.706 
1.294 



0.6159 
0.7972 
0.6563 
0.7678 
0.6664 
0.6644 
0.7398 
0.7995 
0.7111 
0.8068 
0.7087 
0.6739 
0.7041 
0.7556 
0.8645 
0.8734 
0.8822 
0.8911 
0.9000 



0.8120 
0.6442 
0.7601 
0.6736 
0.8694 
0.8380 
0.5724 
0.6931 
0.6387 
0.6189 
0.7067 
0.6626 
0.6394 
0.7901 
0.7286 
0.6510 
0.6100 
0.6421 
0.7313 



0.15 
0.11 
0.16 
0.13 
0.15 
0.16 
0.12 
0.11 
0.12 
0.11 
0.13 
0.12 
0.13 
0.13 
0.12 
0.11 
0.10 
0.09 
0.12 



0.24 
0.18 
0.25 
0.22 
0.23 
0.24 
0.20 
0.18 
0.19 
0.18 
0.21 
0.19 
0.23 
0.19 
0.19 
0.18 
0.17 
0.15 
0.19 



relatively simple, allowing us to work through - in a concrete 
setting - the many steps which will be routinely required in 
the future. In some ways the problem is one of the simplest 
currently confronting theorists, but as we shall discuss below 
even it has demanded collaboration with other communities, 
the development of a significant infrastructure, new modes of 
working, and large amounts of manpower and computational 
capacity. 

In Paper I of this series (Heitmann et al. 2009a) we demon- 
strated that it is possible to obtain nonlinear matter power 
spectra with percent level accuracy out to k ~ 1 h Mpc -1 and 
derived a set of require ments for such simulations. Paper II 
dHeitmann et al. 2009b) described the construction of an em- 
ulation scheme to predict the nonlinear matter power spec- 
trum and the underlying cosmological mod els, building on 
the "Cosmic Calibration Framework" (|Heitmann et al.1 120061: 
lHabib et al.ll2007b ISchneider et al. 2008). Here, in Paper III, 
we present results from the complete simulation suite based 
on the cosmologies presented in the second paper and pub- 
licly release a precision power spectrum emulator 1 . The sim- 
ulation suite is called the "Coyote Universe" after the cluster 
it has been carried out on. We will extend our work to include 
other measurements, such as the mass function or higher order 
statistics, in future publications. 

The outline of this paper is as follows. In ^2] we describe 
the simulations we have run in support of this program, and 
the codes with which they were run. ^describes how we put 
together the different estimations of the power spectra from 
our multi-scale runs while ^presents the details of the emu- 
lator and tests. We describe some lessons learned in §|5]before 
concluding in §|6] 

2. THE SIMULATION SUITE 

2.1. The Simulations 

The Coyote Universe simulation suite encompasses nearly 
1,000 simulations of varying force and mass resolution. The 
simulation volume is the same in all cases, a periodic cube of 
side length 1300 Mpc. We consider 37+1 cosmological mod- 

1 http://www.lanl.gov/projects/cosmology/CosmicEmu 



els, listed in Table Q] which we select with two aspects in 
mind: our statistical framework and current constraints from 
a variety of cosmological measurements (see Paper II for fur- 
ther discussions and see below for a short summary). 

The 37 models, labeled 1-37, are used to construct the em- 
ulator while model M000 is used as an independent check on 
the power spectrum accuracy in the parameter regime of most 
interest. (Other tests are described in §|4]and in Paper II.) 

For each cosmology we run 20 realizations, 16 lower res- 
olution simulations covering the low-A: regime (which we re- 
fer to as the "L-series"), 4 medium resolution runs to provide 
good statistics in the quasi-linear to mildly nonlinear regime 
(the "H-series"), and one high resolution run to extend to 
k~ 1/jMpc" 1 (the "G -series"). The high-resolution run uses 
the same realization as one of the medium resolution runs. 
The low- and medium-resolution runs are performed with a 
particle-mesh (PM) code. The code evolves 512 3 particles in 
the low -resolution runs and 1024 3 particles in the medium- 
resolution runs, in each case with a force mesh twice as large 
in each dimension as the particle load. Densities and forces 
are computed using Cloud-in-Cell (CIC) interpolation and a 
Fast Fourier Transform (FFT)-based Poisson solver. The po- 
tential is determined from the density using a 1/k 2 Green's 
function and the force is computed by 4 th order differencing. 
Particles are advanced with second order (global) symplectic 
time-stepping with Alna = 2%. In order to resolve the high-£ 
part of the power spectrum, we evolve one of the 1024 3 par- 
ticle initial distr ibutions with the Tree-PM code GADGET-2 
(Springel2005). We use a PM grid twice as large, in each di- 
mension, as the number of particles, and a (Gaussian) smooth- 
ing of 1.5 grid cells. The force matching is set to 6 times the 
smoothing scale, the tree opening criterion to 0.5% and the 
softening length to 50kpc. The starting redshift for each sim- 
ulation is z = 211. Further details regarding the simulations 
and choices of simulation parameters can be found in Paper I. 

All the simulations are carried out on Coyote, a large HPC 
Linux cluster consisting of 2580 AMD Opterons running at 
2.6 GHz. The low-resolution PM runs are carried out on 64 
processors, the medium-resolution PM runs and GADGET-2 
runs on 256. For the billion particle runs we store particle and 
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halo information at 1 1 different redshifts between z = 4 and 
Z = 0. This leads to a final simulation database of size roughly 
60TB. 



2.2. The Cosmological Models 

The selection of the cosmological model suite depends on 
two considerations: the statistical framework we use to con- 
struct the emulator and cur rent parameter constra ints from the 
CMB as set by WMAP-5 dKomatsu et al.ll2008h . We do not 
insist on a formal methodology to make the model selection, 
but instead apply some practical and conservative arguments 
to justify our decisions. For an in-depth discussion we re- 
fer the reader to Paper II. In this paper we briefly summarize 
some of the considerations since these will define the region 
for which the emulator will be valid. 

Our aim is to find a distribution of the parameter settings 

- the simulation design - which provides optimal coverage 
of the parameter space, using only a limited number of sam- 
pling points. Simulation designs well-suited for this task are 
Latin-Hypercube (LH)-based designs, a type of stratified sam- 
pling scheme. Latin hypercube sampling generalizes the Latin 
square for two variables, where only one sampling point can 
exist in each row and each column. A Latin hypercube sample 

- in arbitrary dimensions - consists of points stratified in each 
(axis-oriented) projection. 

Very often LH designs are combined with other design 
strategies such as orthogonal array (OA)-based designs or are 
optimized in other ways, e.g., by symmetrizing them (more 
details below). By intelligently melding design strategies, dif- 
ferent attributes of the individual sampling strategies can be 
combined to lead to improved designs, and shortcomings of 
specific designs can be eliminated. As a last step, optimiza- 
tion schemes are often applied to spread out the points evenly 
in a projected space. One such optimization scheme is based 
on minimizing the maximal distance between points in the pa- 
rameter space, which will lead to more even coverage. Two 
design strategies well suited to cosmological applications in 
which the number of parameters is much less than the num- 
ber of simulations that can be performed are optimal OA-LH 
design strategies and optimal symmetric LH design strategies. 
For this project, we have generated forty different designs 
following different strategies and chose the best design from 
these (where "best" refers to best coverage in parameter space 
with respect to specific distance criteria explained in Paper II). 

In order to restrict the number of necessary simulation runs 
to an as small number as possible, it is helpful to keep the 
number of cosmological parameters and their prior ranges 
small. Current observations of the CMB and the large scale 
structure are consistent with a ACDM model with constant 
dark energy equation of state, w. We therefore concentrate 
on the following five cosmological parameters: u) m = Q m h 2 , 
ujf, = ftbh 2 , n s , w, and a% where fl m contains the contributions 
from the dark matter and the baryons. We restrict ourselves 
to power-law models (no running of the spectral index) and 
to spatially flat models without massive neutrinos. For each 
cosmology, h is determined by the angular scale of the acous- 
tic peaks in the CMB (Paper II) which is known to very high 
accuracy (0.3%: iKomatsu et al.ll2008l) . From WMAP 5-year 



data, in combination with BAO, we have 2 

u m = 0.1347 ±0.0040 (3%), 

uj b = 0.0227 ±0.0006 (3%), (1) 

n s = 0.9610 ±0.0140 (2%). 

Current data restrict a constant equation of state for the dark 
energy to w = -1 with roughly 10% accuracy and recent deter- 
minations put the normalization in the range 0.7 < erg < 0.9 
with still rather large uncertainties. 

Considering all these constraints and their uncertainties, we 
choose our sample space boundaries for the 37+1 models to 
lie within the range(s) 

0.120 < u m < 0.155, 
0.0215 < uj b < 0.0235, 

0.85 <n s < 1.05, (2) 
-1.30 <w< -0.70, 
0.61 <cr 8 <0.9, 

over which the emulator is designed to produce reliable re- 
sults. We verified in Paper II that 37 models spanning these 
parameter ranges are indeed enough to generate an emula- 
tor at the 1% accuracy. We emphasize that our emulator is 
valid for the complete parameter space defined by these priors 
and not restricted to some values in a band around the best- 
fit cosmology (for current observational data). Obviously, the 
emulator quality will be slightly worse on the edges of the 
hypercube but the emulator is always accurate within ~ 1% 
anywhere within the priors. 

2.3. Power Spectra 

In Paper I we describe in detail how we obtain the mat- 
ter power spectrum from a snapshot of the simulation. We 
briefly summarize the salient points here. We compute the 
dimensionless power spectrum, 



which is the contribution to the variance of the density pertur- 
bations per \nk. We obtain A 2 by binning the particles onto a 
2048 3 grid using CIC assignment, applying an FFT, correct- 
ing for the charge-assignment window function, and averag- 
ing the result in fine bins in |k| spaced linearly with width 
Ak ~ 0.001 Mpc -1 . As discussed in Paper I, we do not cor- 
rect for particle discreteness as our particle loading is high 
enough to make such corrections unnecessary and there are 
some indications that a simple Poisson shot-noise form is not 
correct. 

3. POWER SPECTRUM DETERMINATION 

When creating a nonlinear power spectrum emulator, we 
prefer that the underlying training data set be smooth. We de- 
scribe here how we construct a smooth power spectrum from 
the 20 realizations of each cosmology and from cosmological 
perturbation theory. 

In each simulation the modes in the initial density field are 
a single realization of a Gaussian random field and this in- 
troduces large run-to-run scatter. This scatter is reduced at 
higher k by the relatively large number of modes which are 
averaged. However, at low k the estimates of the power spec- 
trum exhibit significant scatter which is expensive to reduce 

2 See http://lambda.gsfc.nasa.gov 
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FIG. 1. — Comparison of perturbation theory with the nonlinear matter 
power spectrum from simulations at z = 0. Out of our 38 models we show 
three random examples, the results for the remaining models being very sim- 
ilar. For each model we determine a simulated power spectrum by matching 
the sixteen lower resolution L runs to the four medium resolution H runs at 
k = 0.25 Mpc" 1 . Therefore we average over 20 (1.3Gpc) 3 simulations at low 
^-values and four simulations beyond k = 0.25 Mpc~' . We then take the ratio 
of the simulation results with respect to perturbation theory. The agreement 
is very good out to at least k ~ 0.03 Mpc -1 . We therefore conclude that 
- up to these wavenumbers - perturbation theory results are robust for the 
cosmologies considered in this paper. 

by brute force, i.e. running a very large number of realizations 
in large volumes. In addition, the approach to linear theory at 
low k can be quite slow for many currently popular models 
around ACDM (if 1 % accuracy is the desired goal) so simply 
replacing the N-body results with the input linear model can 
be relatively inaccurate. This is however an area where pertur- 
bation theory can be of help, since the real-space mass power 
spectrum is computable in perturbation theory and there is a 
small but non-negligible range of scales where perturbation 
theory improves upon linear theory. For a recent overview of 
the performance of diffe rent perturbation theory approaches, 
see lCarlson et al.l (I2009I) . 

3.1. Perturbation theory 

While there has been a resurgence of interest in pertur- 
bative methods in recent years, we stick to the simplest 
and oldest method: "standard perturbation theory" ([P eebles 
1980t Uuszkiewicz|[l98Tt lVishniaclll98l iGoroff et all 1 19861; 
Makin oet alj|1992t iJain & Bertschinger 1 1 9941) . We consider 
only the first correction to linear theory, which in standard 
perturbation theory can be written in terms of a simple inte- 
gral over the l inear theory power spectrum (we u se the specific 
form given in Meik sin. White & Peacockll 19991) . When com- 
pared to our simulations, we find that standard perturbation 
theory is accurate at the percent lev el for k < 0.5 k n \ where k n \ 
can be defined as (Matsubara 2008|): 

, 1 f dk A 2 (k) 

^-JT-W- (4) 

The values for k n \ at z = and z = 1 are listed for all models 
in Table [1] For example, for model M000, k „\ ~ 0.1 Mpc" 1 at 
Z = 0. Similar results hav e been reported in Matsubara (120081) 
and lCarlson etaTI J2009I) . 

Almost any scheme that switches smoothly from standard 
perturbation theory to our N-body results around 0.5 £ n i pro- 



duces results that agree at the percent level. At z = 0, this 
would lead to a matching point between k = 0.045 Mpc -1 
(M036) and k = 0.075Mpc _1 (M019). For simplicity we keep 
the matching point the same for all cosmologies. Since we 
have good statistics from our simulations at k ~ 0.03 Mpc -1 
already, we choose this k value as a very conservative match- 
ing point for all models. To verify this point, we compare 
our simulation results to perturbation theory for the different 
cosmologies. The simplest approach for this comparison is 
"brute force": simply take the ratio of the simulation result to 
the perturbation theory prediction. The major obstacle here is 
the run-to-run scatter in the simulations. In order to overcome 
this problem, one can follow two routes: either incorporate 
fluctuations from the realization into the pe rturbation theory 
prediction, as done in Takahas hi et al.l (120081) . or average over 
a large number of large volume simulations. We follow the 
second approach to avoid any ambiguities and systematic er- 
rors resulting from finite box size effects. 

An example of our approach is shown in Figure [T]for a ran- 
dom subset of 3 models from the total of 37 used to build 
the emulator. Although at the lower end of the k range, 
the large but finite sampling volume becomes clearly evi- 
dent, the matching to perturbation theory at k ~ 0.03 Mpc -1 
works extremely we ll. A similar result can also be found in 
ICarlson et al]d2009l) . 

3.2. The Estimation Procedure 

Next, we show how we can combine perturbation the- 
ory and results from A^-body simulations to generate smooth 
power spectra which will be the foundation for building our 
emulator. Two problems have to be solved for this: (1) We 
have to eliminate the scatter in the A^-body results for the 
nonlinear power spectrum without erasing subtle features and 
match results from different resolution simulations. (2) We 
have to match very accurately between perturbation theory 
and simulation results. In the following, we discuss an ap- 
proach based on process convolution to solve these problems. 

3.2.1. Power Spectrum Estimation using Process Convolution 

In this section, we discuss the procedure for estimating the 
smooth power spectrum for each cosmology based on the sim- 
ulation results which possess inherent scatter. Figure|2]shows 
the power spectrum from the simulations for model MOOl at 
z = 0. The data have three notable features that are important 
for the modeling procedure. (1) The non-standard representa- 
tion for the power spectrum 

V(k) = A 2 (k)/k LS (5) 

is chosen to accentuate the baryon acoustic oscillations. We 
will account for this feature when we choose our function 
class for the smooth power spectra (discussed later). (2) The 
three series of simulations; G, H, and L with 1, 4, and 16 
realizations respectively. Because the H and L series do not 
have enough force resolution to resolve the nonlinear regime 
at high k, we need to restrict the use of these runs to small 
and intermediate k ranges. (3) The simulation variance at any 
given k is known. 

We treat each simulated ("noisy") realization of the (large 
volume or averaged) spectrum as a draw from a multivariate 
Gaussian distribution whose mean is given by an unknown 
smooth spectrum. Thus, for a given cosmology c, a given se- 
ries s G {G, H, L}, and a given replicate i = 1 , ■ ■ ■ , N s (where N s 
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FIG. 3. — The process convolution model treats a smooth function as aris- 
ing from a weighted average of a simple stochastic process. In this figure, 
each point on the smooth function is a weighted average of the Gaussian 
impulses shown as the vertical bars. In these case, the weight function is a 
Gaussian kernel. Typically, the smooth function is observed (possibly with 
noise) and the challenge is to estimate the impulses and perhaps some aspect 
of the smoothing kernel. 



FIG. 2. — Power spectra from the W-body simulations for cosmology M001 
on the scales that are used for smoothing. The insufficient force resolution of 
the PM runs (H and L) is apparent at large k. To generate a smooth power 
spectrum over the whole k range, we match perturbation theory and the L- 
runs at k = 0.03 Mpc~', the L- and H-series at k = 0.25 Mpc~', and use the 
G-series for k > 0.32 Mpc~' . Note the non-standard representation of the 
power spectrum via A. 2 (k)/k l 5 . 



is the number of simulations for the series), we have a multi- 
variate Gaussian density for the simulated spectrum P c si , 



f(P c J<x\A s QA] 



1 1/2 



(6) 



x exp 



-A S V C ) 



A S QA T S (P c s j -A s V c 



Here, V c is the smooth power spectrum; A s is a projection 
matrix of zeros and ones that is used to remove the high-A; 
values for which the H and L series are not used (thus, A c is 
an identity matrix); and ft is a diagonal matrix of the known 
precisions (inverse variances). 

The model is completed by specifying a class of functions 
for the smooth power spectrum V°. We choose a flexibl e class 
of sm ooth functions called a process convolution (Higdon 
2002). These functions are best described constructively. A 
process convolution builds a smooth function as a moving 
average of a simple stochastic process like independent and 
identically-distributed (i.i.d.) Gaussian variates or Brownian 
motion. The moving average uses a smoothing kernel whose 
width is allowed to vary over the domain to account for non- 
stationarity (smoother in some regions, more wiggly in oth- 
ers). Figure [3] shows a simple example of a process convolu- 
tion built on white noise smoothed with a Gaussian kernel. 

We build the process convolution for V c on Brownian mo- 
tion u c , with marginal variance t%, realized on a sparse grid 
(relative to the power spectrum) of evenly spaced points, x 
(the number of points is not important so long as it is large 
enough; we use 100), 



7 W 



1/2 



exp 



u c 'Wu c 



(7) 



where W is the Brownian precision matrix with diagonal 
equal to [1 2 • • • 2 1] and -1 on the first off-diagonals (this 
matrix cannot actually be inverted, but never has to be in the 
estimation). The Brownian motion is transformed into V c by 



the smoothing matrix K a , 



V c = K a u c 



(8) 



The smoothing matrix K a is built using Gaussian kernels 
whose width varies smoothly across the domain. Thus, we 
have 

i rjogKgo^ (9) 



2a 



where k, is the z'-th value of k for which the power spectrum is 
computed. 

In the description of K a , a is indexed to indicate that it 
changes over the domain. Intuitively, we want a to be small 
in the middle of the domain in order to capture the oscilla- 
tions, but large elsewhere to smooth away the noise. In order 
to estimate this varying bandwidth parameter, we build a sec- 
ond process convolution model. This model is built on i.i.d. 
Gaussian variates v, with mean zero and variance r t 2 , observed 
on an even sparser grid of evenly spaced points, t (length M v ; 
we use 10, but any large enough number will suffice), 



/(v)« 



1 



exp 



2r,2 



(10) 



The process v is transformed into a by the smoothing matrix 

K s , 



<7 = K°V. 



(11) 



This matrix is also built using Gaussian smoothing kernels, 
but with a constant bandwidth, S. Thus, we have 



Kf; = 



exp ■ 



2S 2 



(12) 



Combining all of this, we get a distribution for the simu- 
lated power spectra for a given cosmology, 



,/,l/2 



(13) 



/(/£)cx|A s fiA; 

x exp j-I (F Sii -A s K°u c )'A s £lA' s (i^-A.W) J , 

where s e {G, H, L} and i=l,-~,N s . Note that K a is a func- 
tion of the parameters v and 5. We choose noninformative 
priors for t%, t, 2 , and 5 so as to impart little or no information 
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FIG. 4. — Predictions for the linear power spectrum from 20 realizations of the initial power spectrum for the same three models shown in FigurefT] The upper 
panels for each model show in gray the realizations, in red the smooth power spectrum and in black (dashed) the linear theory power spectrum that underlies the 
simulations. The lower panels show the realizations (again in gray) divided by the linear theory answer and in red the smooth power spectrum divided by linear 
theory. The blue line shows the empirical mean of the realizations divided by linear theory. The vertical red lines in each panel indicate the matching point of 
linear theory and simulation result, the horizontal dashed lines in the lower panels show the 1 % error. In all models the discrepancy is at the 1 % level. 



about their values, 



7r(T„ ) OC 



7r(T v ) (X 



exp 



exp 




(14) 



n(S)(xI{0<d < 10}, 



which are inverse gamma, inverse gamma, and uniform, re- 
spectively. 

Equations (0 (for all c), (fTDt. and ([141 are multiplied 
together with Eqn. dT3l> (for all c) to produce a posterior 
distribution for the unknown parameters and stochastic pro- 
cesses. Obtaining an estimate for the smooth power spec- 
tra requires an estimate for each of the parameters, the pro- 
cess v, and each of the processes u c for all c. Markov Chain 
Monte Carlo (MCMC) vi a the Metroplis-Hastings algorithm 
dChib & Gr eenberg 1995) produces a sample from the poste- 
rior distribution by drawing each parameter individually. The 
stochastic processes u c can be integrated out of the distribu- 
tion, so the MCMC produces samples of the parameters as 
well as v. We use the posterior mean of the parameters to 
obtain a conditional mean for each of the u c which is then 
transformed into each of the V°. 

The perturbation results are included by setting the low k 
values for each of the simulated spectra in a given cosmology 
to the perturbation results and setting the precisions for these 
values to be very large. This is done prior to the estimation 
described here. The replication of these values in every sim- 
ulation, combined with the large precisions, nearly forces the 
estimated result through these points. Further, the transition 
from the N-body results to the perturbation results will be rel- 
atively smooth as long as the two line up fairly well. 

3.2.2. Tests with the Linear Power Spectrum 

The next step is to test the matching procedure between the- 
oretical and simulation results to ensure high accuracy at the 
perturbation theory matching point. Since we do not know the 
exact answer for the nonlinear power spectrum, we first carry 
out a test using linear theory. In this test we use the power 
spectra from the initial conditions and show how well we can 
smooth out the run-to-run scatter. Knowing the exact answer 



here will allow us to assess how well the matching procedure 
actually works. 

For these spectra, the entire estimation procedure is applied 
to the initial condition spectra. The only difference is that 
fewer grid points were used for the latent processes u and v to 
account for the reduced range of the spectra (we use 70 and 7, 
respectively). The final results are shown in Figure |4]for the 
same set of three random models as in Figure Q] We verified 
that the results hold for all of the remaining models. The up- 
per panels of each sub-panel show the theoretical linear power 
spectrum in black and the prediction from the simulations in 
red. The vertical line marks the matching point to linear the- 
ory at k ~ 0.03Mpc _1 . The lower panels show the ratio of 
the predicted power spectra to the theoretical power spectra 
in red. Below the matching point, the agreement is - by con- 
struction - perfect. Beyond the matching point, the smoothed 
prediction in red is accurate at the 1 % level. This test shows 
that we can obtain a smooth, high-accuracy prediction for the 
power spectrum by combining a suite of realizations with per- 
turbation theory. 

3.3. The Nonlinear Power Spectrum 

The process convolution procedure is applied to the 
power spectrum realizations for each of the 37 cosmolo- 
gies at six values of the scale factor a = 1/(1+ z) € 
{0.5,0.6,0.7,0.8,0.9,1.0}, where z is the redshift. The re- 
sults are shown in Figure|5] again for a subset of three models. 
Although we have no known truth to use as a comparison in 
this case, the resulting estimates continue to fit the simulation 
realizations very well. 

In place of comparing with known truth, we can examine 
some tests of our modeling assumptions. We assume that the 
simulation spectra on the modified scale are independently 
and normally distributed about a smooth mean with a variance 
that changes with k. Given this, we can compute standardized 
residuals in which the estimated smooth mean is subtracted 
from the simulations and the result is scaled by the known 
standard deviation. These standardized residuals should look 
like i.i.d. standard normal variables. Figure|6]shows quantile- 
quantile plots for the G simulations of three cosmologies at 
each of the six scale factors. Theoretical quantiles from a 
standard normal distribution are given on the jt-axis and the 
sample quantiles of the standardized residuals are shown on 
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FIG. 5. — Simulated power spectra and the smooth estimate for each of the three cosmologies at six values for the expansion factor a. 



the y-axis. The nearly straight line at 45° indicates an ex- 
tremely good distributional fit. Figure [7] shows the standard- 
ized residuals plotted against k for the same simulations. This 
plot verifies the independence assumption. Further, there is 
an almost complete lack of evidence of any structure in these 
plots, suggesting that we would not greatly improve the fit by 
relaxing the smoothness assumption. 

It is also interesting to examine the plot of the kernel width 
function as estimated by the MCMC process. Figure [8] shows 
the median draw for a as a function of k. As expected, the 
kernel width is small in the vicinity of the baryon wiggles. 
This means the values of the latent process u in this region 
receive large weights, with the contributions for values further 
away dropping off very quickly. On both ends, the kernel 
width is large, which results in a smooth function in these 
regions. 

4. THE EMULATOR 

Having extracted the smooth power spectra from our sim- 
ulation suite, we can now build an emulator to predict the 
nonlinear matter power spectrum within the priors specified 
in Eqns. (|2). We will use only models 1-37 for the emulator 
construction; model M000 will serve as an independent check 
of the emulator accuracy, along with hold-out tests described 
below. 

In order to construct the emulator, we model the 37 power 
spectra using an n-p dimensional basis representation: 

Up 

V(k;z;6) = J2<f> i (k;z)w i (0), 6e[0,l]"\ (15) 

where the <fii(k;z) are the basis functions, the w,-(0) are the 
corresponding weights, and the 9 represent the cosmological 
parameters. The dimensionality np refers to the number of 
orthogonal basis vectors {</>,(£, z), . . . : 4> nv {k 1 z)}- The param- 
eter rig is the dimensionality of our parameter space - with 5 
cosmological parameters we have rig = 5. The power spectrum 
V(k\ z; 0) depends on the wavenumber k, the redshift z, and the 
five cosmological input parameters 8 (note that we rescale the 
range of each parameter to [0, 1]). Examination of the results 
indicates that n-p = 5 is a good choice for the number of basis 
vectors </>,(£; z) (that np = tig here is a coincidence). The task is 
now to (1) construct a suitable set of orthogonal basis vectors 
4>i(k;z) and (2) model the weights Vfj(fJ)). For the first task we 



use principal components, and for the second, Gaussian Pro- 
cess (GP) models. Our choice of GP modeling is based on 
their success in representing functions that change smoothly 
with parameter variation, e.g., the variation of the power spec- 
trum as a function of cosmological parameters. Both steps are 
explained in great detail in Paper II; we refer the interested 
reader to that publication. Paper II also contains various error 
control tests of the GP-based interpolation method which we 
do not repeat here. 

Since the details of how to build an emulator are already 
provided in Paper II we can immediately turn to our final 
product: the emulator itself. To facilitate use of the emulator, 
we are releasing the fully trained emulator with this paper. It 
provides nonlinear power spectra at a set of redshifts between 
z = and z = 1 out to k = 1 /iMpc" 1 , for any cosmological 
model specified within the priors given in Eqns. (0. 

4.1. Emulator Test 

In order to verify the accuracy of the emulator, we perform 
two important checks (other tests on the methodology were 
carried out in Paper II). For the first, we compare the simu- 
lation results for model M000 with the emulator prediction. 
Figure [9] shows the results of this true out-of-sample predic- 
tion. The plot shows the ratio of the emulated power spectrum 
to the actual simulation for the M000 cosmology at six values 
of the scale factor a. This cosmology is completely interior 
in the design and, since M000 was not used to estimate the 
smoothing or the emulator, this test provides a good indication 
of the actual performance of the emulator within its parame- 
ter bounds. Overall the agreement between the simulations 
and the emulator is excellent and for most of the fc-range well 
below the 1 % error bound. 

The second check consists of a sequence of holdout tests. 
In a holdout test, the emulator is built from 36 cosmological 
models and the emulator prediction can be then compared to 
the result from the extra model. The drawback of this test is 
that if we have only a very small number of models, each of 
them is important for capturing some part of the parameter 
space and taking it out for building the emulator degrades the 
emulator precision. In order to keep this problem to a min- 
imum, we only perform holdout tests for models which are 
interior simulations, meaning that none of the five parameters 
are near the extreme limits of the chosen prior range. There 
are six such models: M004, M008, MOB, M016, M020, and 
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FIG. 6. — Quantile-quantile plots of the standardized residuals for the G simulation at the six scale factors for three cosmologies. Standardized residuals are 
computed by subtracting the estimated mean from the simulation and multiplying each value by the square root of the known precision at its k value. Our 
assumptions suggest that the resulting sample should follow a standard normal distribution with no dependence on k. These plots show the sample quantiles of 
the standardized residuals (essentially the sorted values) plotted against the theoretical quantiles for a sample of this size from the standard normal distribution. 
The nearly straight lines indicate little deviation from our assumptions and suggest that the model fits well. 
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FIG. 7. — The standardized residuals for the G simulation at the six scale factors for three cosmologies plotted against k. There are no obvious correlations or 
structure confirming our assumptions. 



M026. The ratio of the emulated prediction to the actual sim- 
ulation for these models is shown in Figure [10] For each cos- 
mology, there are six ratios for each value of the scale factor. 
The lines are quite well behaved with errors largely within the 
1% bounds. 

5. LESSONS LEARNED AND FUTURE CHALLENGES 

The advent of precision cosmology and the prospect of very 
large surveys such as LSST and JDEM poses an enormous 
challenge to the theory community in the field of large scale 
structure predictions. Future progress not only requires very 
accurate predictions but also the ability to produce predic- 
tions for different cosmological models very fast. The aim of 
the Coyote Universe project was to take a first step in attack- 
ing this problem, focusing on the matter power spectrum at 
intermediate scales out to k ~ 1 /zMpc" 1 . While predicting the 
matter power spectrum on these scales at high accuracy may 
appear to be a moderately difficult task, it turned out to be a 
technical and computational challenge in several respects. It 
is generally hard to imagine problems and pitfalls in advance, 
which is why fully working through an example is so help- 
ful. Since the community has to follow a similar path in the 



future to create predictive capabilities for different cosmolog- 
ical probes, we summarize here some of the lessons learned 
during this project. 

The major differences between a project like this (including 
all three papers of the series) and previous numerical studies 
of large scale structure probes are: 

(1) Computational and Storage Capacity: The Coyote Uni- 
verse simulation suite encompasses roughly 60 TB of data. 
The computational cost is of the order of a million CPU-hours 
and including waiting times in submission queues, downtimes 
of the machine, and so on, carrying out these simulations took 
roughly six months. The simulation size of one billion parti- 
cles in a Gpc 3 volume was barely enough to resolve the scales 
of interest and for future work will certainly not be sufficient. 
Such simulations will need larger volumes and better mass 
resolution. For example, to resolve a 1O 12 M halo with 100 
particles in a (3 Gpc) 3 volume, we would need 300 billion par- 
ticles. While supercomputers will get faster and larger in the 
future, generating many simulations at the edge of machine 
capabilities will always be a challenge. In addition, archiving 
the outputs of such simulations will become very expensive 
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FIG. 8. — Median MCMC draw for the bandwidth function a (black line). 
Small values on this plot correspond to places where the spectra are compar- 
atively less smooth because of the baryon wiggles. In addition, all the power 
spectra are shown (37 models, 6 redshifts each) in light gray. The power 
spectra a scaled and shifted to fit the plot, cr is low on the baryon acoustic 
oscillation scale to accommodate local wiggles, as expected. The vertical 
line shows the matching point between perturbation theory and the simula- 
tion outputs. On the left of this, each replicate is identical (compared to the 
realizations from the different simulation boxes) and smooth, therefore the 
behavior of a has not much information. 
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FIG. 9. — Ratio of the emulator prediction to the smooth simulated power 
spectra for the M000 cosmology at six values of the scale factor a. The error 
exceeds 1 % very slightly in only one part of the domain for the scale factors 
« = 0.7,0.8, and 0.9. 

in terms of storage. From the Coyote Universe runs we stored 
11 time snapshots plus the initial conditions (particle posi- 
tions and velocities and halo information) leading to 250 GB 
of data per run. For the 300 billion particle run this would 
increase to 75 TB. Only very few places worldwide would be 
able to manage the resulting large databases. 

(2) Simulation Infrastructure: Running a very large num- 
ber of simulations makes it necessary to integrate the ma- 
jor parts of the analysis steps into the simulation code and 
to automate as much of the mechanics of running the code 
(submission, restarts) as possible. For the Coyote Universe 
project we developed several scripts to generate the input files 




k[1/Mpc] 

FIG. 10.— M004, M008, MOB, M016, M020, and M026. For each sim- 
ulation, we show six residuals corresponding to the different values of the 
scale factor a. The errors are on the order of 1% for the bulk of the domain 
of interest. Considering that the tested emulators are built on an incomplete 
design, this result is remarkably good. 

of the codes, to structure the directories where different runs 
are performed in, and to submit the simulations to the comput- 
ing queue system. For future efforts of this kind the adoption 
and development of dedicated workflow capabilities for these 
tasks must be considered. The number of tasks to be car- 
ried out will become too large to keep track of without such 
tools. In addition, since large projects will require extensive 
collaborations, software tools will make it easier to work in a 
team environment since each collaborator will have informa- 
tion about previous tasks and results. An example of such a 
tool for cosmological simulatio n analysis and visualization is 
given in lAnderson et alJ (120081) . 

We carried out the data analysis after the runs were finished. 
For very large simulations this is not very practical, and on- 
the-fly analysis tools are required to minimize read and write 
times and failures. This in turn requires that the code infras- 
tructure be tailored to the problem under consideration. 

(3) Serving the Data: Clearly, large simulation efforts can- 
not be carried out by a few individuals, and require possi- 
bly community wide coordination. The simulation data will 
be valuable for many different projects. It is therefore nec- 
essary to make the data from such simulation efforts pub- 
licly available and serve them in a way that new science can 
be extracted from different groups of simulations. Transfer- 
ring large amounts of data is difficult because of limitations 
in communication bandwidth and also because of the large 
storage requirements. It would therefore be desirable to have 
computational resources dedicated to the database. In such 
a situation, researchers would be able to run their analysis 
codes on machines with direct access to the database and per- 
form queries on the data easily. We are planning to make the 
Coyote Universe database available in the future and use it as 
a manageable testbed for such services. 

(4) Communication with other Communities: The complex- 
ity of the analysis task makes it necessary to efficiently collab- 
orate and communicate with other communities, for example 
statisticians, computer scientists, and applied mathematicians. 
Many tools that will be essential for precision cosmology in 
the future have already been invented - the task is to find them 
and use them in the best way possible. 

6. CONCLUSIONS 
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This paper is the last of the Coyote Universe series of pub- 
lications. Paper I was concerned with demonstrating that per- 
cent level accuracy in the (gravity only) nonlinear power spec- 
trum could be attained out to k ~ 1 Mpc -1 . Paper II showed 
that with only a relatively small number of simulations, in- 
terpolation across a high-dimensional space was possible at 
close to the same level of accuracy as that attained in the in- 
dividual runs. Paper III takes this work to the final conclu- 
sion: Based on almost 1000 simulations spanning 38 wCDM 
cosmologies, we present a fast and very accurate prediction 
scheme - an emulator - for the nonlinear matter power spec- 
trum. The emulator is accurate at the percent level, improving 
over commonly used fitting functions by almost an order of 
magnitude. 

The emulator construction - as explained in Paper II - is 
based on GP modeling. In order to carry this out, a major 
challenge is to produce a smooth power spectrum from a finite 
set of simulations for each cosmology. In order to minimize 
run-to-run scatter on very large scales, we performed several 
medium resolution simulations and matched these at suffi- 
ciently low k to perturbation theory. On quasi-linear scales 
we used medium resolution simulations and matched those to 
high resolution runs at small spatial scales. Matching the dif- 
ferent resolution runs and perturbation theory accurately was 
carried out using process convolution. This technique allowed 
us to construct a smooth power spectrum for each cosmologi- 
cal model. The results were then used to construct the emula- 
tor via GP modeling as described in detail in Paper II. 

We are releasing the power spectrum emulator as a C code, 
which allows the user to specify a cosmology within our pri- 
ors and returns the power spectrum at six different redshifts 
between z = and z — 1 out to k ~ 1 Mpc" 1 . These power 
spectra can now be used for further analysis of cosmological 
data. We are planning to extend the emulator in the near future 
to a larger Arrange and will provide a smooth interpolation be- 
tween results from different redshifts. 

A major challenge will be to ensure that a certain level 
of accuracy is reached with the simulations when going to 
small scales. Besides being computationally very expensive 
(high force and mass resolution being required) the physics 
at smaller scales is far more complicated. The inclusion of 



gasdynamics and feedback effects (along with other physics) 
is far from being straightforward. The impossibility of car- 
rying out a direct simulation effort is more or less certain; a 
good number of phenomenological/subgrid modeling param- 
eters will be required. As simulation complexity and the num- 
ber of modeling and cosmological parameters increases, it be- 
comes even more important to develop efficient and controlled 
sampling schemes as described in Paper II, so that data can be 
used to determine both cosmological and modeling parame- 
ters (self-calibration). 

With the series of three Coyote Universe papers we have 
demonstrated that it is possible to extract cosmological statis- 
tics such as the power spectrum at high accuracy and that one 
can build an accurate prediction scheme based on a limited 
set of simulations. This line of work will be important for in- 
terpreting results of future cosmological surveys. It will also 
have to be extended in several ways: (1) the cosmological 
model space has to be opened up; (2) we have to ensure high 
accuracy at scales smaller than those considered here; (3) we 
have to include more physics in order to capture those small 
scales correctly; (4) we have to include different cosmologi- 
cal probes, e.g., the cluster mass function and the shear power 
spectrum, to be able to build a complete framework for an- 
alyzing future survey data. We have shown here that such a 
program can in principle be established, though it will demand 
a large concerted effort between different communities. 
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